program main
use type_kinds
implicit none

double precision::omega,delta_t,delta,Rsave,Rin1,Rin2,Rin3,R,T,eps,s
double precision, dimension(:,:),allocatable::h,update,transp,h_exist
integer::i,j,k,numrow,numcol,iterations,p_count,pump_on,out_count,day_count

open(unit=105,file='aqout_days.txt',status='replace')
open(unit=103,file='aqout_all.txt',status='replace')
open(unit=104,file='aqout_select.txt',status='replace')
open(unit=200,file='aqout_downdraw.txt',status='replace')

!==============================================================!
Rin1=-1500.0d0						!Pumping Rate (700,400) pump 1
Rin2=-1000.0d0                      !Pumping Rate (2400,2100) pump 2
Rin3=-2000.0d0                      !Pumping Rate (2400,2100) pump 2
omega=1.0
T=200.0d0							!Transmissivity
eps=0.000010d0						!epsilon
s=0.0050
numrow=31							!number of rows
numcol=31							!number of columns   
!==============================================================!

!Start violation days output file
!==============================================================!
write(200,'(T1,a30)')"drawdown (meters),day,position"
!==============================================================!

allocate(h(numrow,numcol))
allocate(h_exist(numrow,numcol))
allocate(update(numrow,numcol))
iterations=0
h=1d0
update=1d0
h(1,:)=50.0d0
h(31,:)=45.0d0

!first calculate the heads of the existing condition.
!==============================================================!
Rin1=0.0d0
Rin2=0.0d0
Rin3=0.0d0
R=0.0d0

!save existing conditions matrix
h_exist=h

do
!always set boundary conditions
update(1,:)=50.0d0
update(31,:)=45.0d0

     do i=2,numrow-1

          do j=1,numcol

             if(j==1)then
                update(i,j) = (omega/4d0)*(h(i-1,j)+2d0*h(i,j+1)+h(i+1,j)+(delta**2d0*R/T)-(4d0*h(i,j)))
                h(i,j)=h(i,j)+update(i,j)
             elseif(j==numcol)then
                update(i,j) = (omega/4d0)*(h(i-1,j)+2d0*h(i,j-1)+h(i+1,j)+(delta**2d0*R/T)-(4d0*h(i,j)))
                h(i,j)=h(i,j)+update(i,j)
             else
                update(i,j) = (omega/4d0)*(h(i-1,j)+h(i,j-1)+h(i,j+1)+h(i+1,j)+(delta**2d0*R/T)-(4d0*h(i,j)))
                h(i,j)=h(i,j)+update(i,j)
             endif
          enddo

     enddo

            if(maxval(abs(update(2:30,:)))<eps)exit

enddo
  
do i=1,numrow
	do j=1,numcol-1
		write(102,'(f6.2)',advance='no')h(i,j)
	end do
	write(102,'(f6.2)')h(i,numcol)
end do


do i=1,numrow
	do j=1,numcol-1
		write(103,'(f6.2)',advance='no')h(i,j)
	end do
	write(103,'(f6.2)')h(i,numcol)
end do

do i=1,numrow
	do j=1,numcol-1
		write(104,'(f6.2)',advance='no')h(i,j)
	end do
	write(104,'(f6.2)')h(i,numcol)
end do
write(105,'(i4)')iterations

!now calculate for days
delta=10.0
Rin1=-1500.0d0 /200.0d0
Rin2=-1000.0d0 /200.0d0
Rin3=-4000.0d0 /200.0d0

!setup a counter for the 20 day pumping rotation
p_count=1
pump_on=1
out_count=1

!Liebmanns method
!==============================================================!
iterations=1
do
!always set boundary conditions
update(1,:)=50.0d0
update(31,:)=45.0d0

     do i=2,numrow-1

          do j=1,numcol
            if(i==5.and.j==8)then
	      R=Rin1
            elseif(i==21.and.j==25)then
	      R=Rin2
            elseif(i==16.and.j==16)then
	      R=Rin1
	      if(pump_on==0)then
		R=0.0d0
	      end if
	    else
              R=0.0
            endif

             if(j==1)then
                update(i,j) = (omega/4d0)*(h(i-1,j)+2d0*h(i,j+1)+h(i+1,j)+(delta**2d0*R/T)-(4d0*h(i,j)))
                h(i,j)=h(i,j)+update(i,j)
             elseif(j==numcol)then
                update(i,j) = (omega/4d0)*(h(i-1,j)+2d0*h(i,j-1)+h(i+1,j)+(delta**2d0*R/T)-(4d0*h(i,j)))
                h(i,j)=h(i,j)+update(i,j)
             else
                update(i,j) = (omega/4d0)*(h(i-1,j)+h(i,j-1)+h(i,j+1)+h(i+1,j)+(delta**2d0*R/T)-(4d0*h(i,j)))
                h(i,j)=h(i,j)+update(i,j)
             endif
          enddo

     enddo
	    p_count=p_count+1


!sets the pump count days and turns pumping on and off
	    if(p_count == 20)then
		if(pump_on==1)then
		  pump_on=0
		else
		  pump_on=1
		end if
		p_count=0
	     end if
		 
!write out any days where drawdown is greater than 2 meters
if(abs(h(16,16)-h(14,16)).gt. 0.50d0)then
write(200,'(T1,f6.2,t7,a1,t8,i4,t12,a7)')abs(h(16,16)-h(14,16)),",",iterations,", north"
else if(abs(h(16,16)-h(18,16)).gt. 0.50d0)then
write(200,'(T1,f6.2,t7,a1,t8,i4,t12,a7)')abs(h(16,16)-h(18,16)),",",iterations,", south"
else if(abs(h(16,16)-h(16,14)).gt. 2.0d0)then
write(200,'(T1,f6.2,t7,a1,t8,i4,t12,a6)')abs(h(16,16)-h(18,16)),",",iterations,", west"
else if(abs(h(16,16)-h(16,18)).gt. 2.0d0)then
write(200,'(T1,f6.2,t7,a1,t8,i4,t12,a6)')abs(h(16,16)-h(18,16)),",",iterations,", east"
end if

call write_output()


	    if(iterations == 1000)exit
		iterations=iterations+1
end do

stop
contains

subroutine write_output()

if(iterations.le.500)then
  if(out_count == 2 .or. out_count ==4)then
    day_count=1
    out_count=out_count+1
  else
    day_count=0
    out_count=out_count+1
  end if
else if(iterations.le.1000)then
  if(out_count ==4)then
    day_count=1
    out_count=out_count+1
  else
    day_count=0
    out_count=out_count+1
  end if
end if
if(out_count==5)out_count=1

!first always write to the all cases file
do i=1,numrow
	do j=1,numcol-1
		write(103,'(f6.2)',advance='no')h(i,j)
	end do
	write(103,'(f6.2)')h(i,numcol)
end do
!write all under 100 to select output file and days
if (iterations.le.100)then
do i=1,numrow
	do j=1,numcol-1
		write(104,'(f6.2)',advance='no')h(i,j)
	end do
	write(104,'(f6.2)')h(i,numcol)
end do
write(105,'(i4)')iterations

!two cases for 100 thru 500 only write even iterations to select
else if(iterations.le.500 .and. day_count==1)then
do i=1,numrow
	do j=1,numcol-1
		write(104,'(f6.2)',advance='no')h(i,j)
	end do
	write(104,'(f6.2)')h(i,numcol)
end do
write(105,'(i4)')iterations

!if iterations between 500-1000 then write every four iterations
else if(day_count==1)then
do i=1,numrow
	do j=1,numcol-1
		write(104,'(f6.2)',advance='no')h(i,j)
	end do
	write(104,'(f6.2)')h(i,numcol)
end do
write(105,'(i4)')iterations
end if

end subroutine write_output
end program
